#ifndef COMMON_RANDOM_H_
#define COMMON_RANDOM_H_

#include <cstdint>

namespace fermat {
	
class Random {
 private:
  uint32_t _seed;
 public:
  explicit Random(uint32_t s) : _seed(s & 0x7fffffffu) { }
  uint32_t next() {
    static const uint32_t M = 2147483647L;   // 2^31-1
    static const uint64_t A = 16807;  // bits 14, 8, 7, 5, 2, 1, 0
    // We are computing
    //       seed_ = (seed_ * A) % M,    where M = 2^31-1
    //
    // seed_ must not be zero or M, or else all subsequent computed values
    // will be zero or M respectively.  For all other values, seed_ will end
    // up cycling through every number in [1,M-1]
    uint64_t product = _seed * A;

    // Compute (product % M) using the fact that ((x << 31) % M) == x.
    _seed = static_cast<uint32_t>((product >> 31) + (product & M));
    // The first reduction may overflow by 1 bit, so we may need to
    // repeat.  mod == M is not possible; using > allows the faster
    // sign-bit-based test.
    if (_seed > M) {
      _seed -= M;
    }
    return _seed;
  }
  // Returns a uniformly distributed value in the range [0..n-1]
  // REQUIRES: n > 0
  uint32_t uniform(int n) { return next() % n; }

  // Randomly returns true ~"1/n" of the time, and false otherwise.
  // REQUIRES: n > 0
  bool oneIn(int n) { return (next() % n) == 0; }

  // Skewed: pick "base" uniformly from range [0,max_log] and then
  // return "base" random bits.  The effect is to pick a number in the
  // range [0,2^max_log-1] with exponential bias towards smaller numbers.
  uint32_t Skewed(int max_log) {
    return uniform(1 << uniform(max_log + 1));
  }
};
}

#endif
